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Abstract 

This paper is concerned with various problems involving the interplay of asymptotics 
and numerics in the analysis of wave propagation in dissipative systems. A general ap- 
proach to the asymptotic analysis of linear, dissipative waves is developed. We apply 
it to the derivation of asymptotic boundary conditions for numerical solutions on un- 
bounded domains. Applications include the Navier-Stokes equations. Multidimensional 
traveling wave solutions to reaction-diffusion equations are also considered. We present 
a preliminary numerical investigation of a thermo-difTusive model of flame propagation 
in a channel with heat loss at the walls. 


1 Introduction 

The mathematical theory of wave propagation in dissipative systems is rich both in appli- 
cations and challenges. Two important classes of problems are those involving convection 
dominated but viscous fluid flows and those involving traveling wave solutions of nonlinear 
equations of parabolic or incompletely parabolic type. (The two classes overlap, of course, 
in the study of viscous shock profiles.) In this paper we shall discuss a variety of problems 
involving the asymptotic and numerical analysis of dissipative waves with an emphasis on 
the connection between the two. 

The construction of numerical radiation boundary conditions at artificial boundaries 
is an issue of great importance in many branches of fluid mechanics including acoustics, 
aerodynamics, meteorology and oceanography. It is also a prime example of an area where 
asymptotic analysis can contribute to the design of numerical methods. Our approach 
is to study the far- field asymptotics of propagating waves and use it to design accurate 
boundary operators. In Section 2 we develop a general asymptotic analysis of the long range 
propagation of waves in linear, dissipative systems. This is based on the identification of 
dominant wave groups associated with minimal decay rates. In Section 3 we show how to 
use these to construct asymptotic boundary conditions and also give some consideration to 
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the small dissipation limit. An early version of these ideas was applied to an advection- 
diffusion equation in [11]. Error estimates are given in [13] while their application to the 
incompressible Navier-Stokes equations is discussed in [12]. 

Traveling wave solutions to reaction-diffusion equations in one space dimension have 
been extensively studied. In Section 4 we consider a natural multidimensional generalization 
of these: traveling wave solutions in cylindrical domains. Far less is known about such 
solutions. It is particularly difficult to generalize the existence theory as it typically makes 
use of a phase space analysis. We consider the numerical computation of such waves. As 
an example we take a thermo-diffusive model of flame propagation in a channel with heat 
loss at the walls. Some preliminary numerical results are presented. All the work in this 
section has been carried out jointly with Steve Buonincontri of Lawrence Livermore National 
Laboratories. 


2 Asymptotic Expansions for Dissipative Waves 


We consider the following system of partially separable second order differential equations 
in an n-dimensional cylindrical spatial domain, ( x,y ) £ [0,oo) x fi: 


du du ^-1 du 


+ = V £ J-Diip-- 


( 1 ) 


(Here we identify x with yo.) These are supplemented by boundary and initial conditions 
defining a signalling problem: 

«(*,y>0) = 0, (2) 

E o u(0, y,t) = g(y,t), (3) 

du 

To— + 2^ = 0, yedtt, (4) 

where ^ denotes the conormal derivative. For now we simply assume that the m x m 
matrices M , A. Bj> C, D tJ , Eq and T t are functions only of the cross-sectional coordinate y 
and that (1-4) is well-posed. A formal representation of the solution may be obtained using 
Laplace transformations in x and t. This leads to the eigenvalue problem in f h 


d 


v ^w, D 


dv ’i — l dv 

v— + E p y 

i = i 


dyj 


dyj 


+ Qv = 0, 


( 5 ) 


Q — sM + A A + C — V(X Dqo + A V'' ), 

U dy ’ 

(6) 

p j = B } - \V(Doj + Djo). 

( 7 ) 

An expression for u then is: 


u(x,y,t) = M x ,y, 0. 

( 8 ) 

Id A' 


w/(*,y,0= f c l (p)qi{x,y,t - p)dp, 
Jo 

where 

qi(x.y.t) - f r ,i+XMj: vi(y,.s)(ls. 

Znl JC 

( 9 ) 

( 10 ) 


(Here / 6 N if 3i(A/) < 0, J?(s) — <• oc.) 
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2.1 Asymptotic Analysis 

We now use the method of steepest descent to compute asymptotic expansions of qi, valid 
for x large. Evaluating (10) along rays t = 71 , x > 1, we seek points s* such that: 

A {(0 = - 7 , (11) 

8(7) > 0, ( 12 ) 

9(7) = 0. ( 13 ) 

Assuming that for 0 < 7 m <„ < 7 < Imax < 00 a solution to (11) exists we have: 

qi (x,y,t) ~ t,t(y; * ,( « )) - ee Mx,y,t), (14) 

7mtn ^ — — 7max* (^) 

Substituting these into (9) formally yields an approximation of ui for t > 7 min x \ 

i n{x,y y t)~ f 7m,n c((p)<M*>sM ^ p) d P- ( 16 ) 

ax % ) 

In (16) ui is represented as a superposition of wave groups propagating at their group 
velocity. Further simplifications may be obtained by noting that, for problems with dissi- 
pation, each wave group decays exponentially with the rate: 

- K( 7 5* + A ,(«*)) / °- (17) 

For general signal data we expect that the large x behavior will be dominated by the 
wave group with least decay, so we seek to minimize the decay rate as a function of 7 . 
Differentiating with respect to 7 yields: 

K(s‘ + ^(7 + A;)) = 0 , (18) 

d 7 


which by ( 11 ) reduces to: 

»(s*( 7 )) = 0. ( 19 ) 

Assuming (11) defines a curve in s space, we have shown that critical points of the decay 
rate occur as the curve crosses the imaginary axis. This may occur for s — 0 or for complex 
conjugate pairs s* = ±Ki. In the spirit of Laplace’s method a crossing at 7 = 7 / will lead to 
a contribution to the asymptotic approximation of u t computed by replacing <pi in (16) by: 


Here we have: 


e x ° r F(x + -JIT’* ~ P) v °(v)’ 


F(z,t) = e 




( a , 1 ) 3 





A? = 7 is* + A/(s* ), 


( 20 ) 

( 21 ) 


( 22 ) 
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( 23 ) 


\ 1 d\[ f - V 

A ' = 77 {3) = ~ v ' 
\2 / *\ 


(24) 


where s* is evaluated at 7 /. For s* = 0 this is indeed Laplace’s method, and the approxi- 
mation will give the leading order asymptotics. For imaginary s* this is not true pointwise 
due to the oscillatory nature of <fi(. The results of [13] imply, however, that it provides an 
asymptotic approximation in an appropriate norm. 


2.2 Sufficient Conditions for ^*(7/) = 0 

Although our analysis allows for s* ^ 0, we are then faced with the potentially difficult 
problem of locating the points on the imaginary s-axis where the group velocity is real. 
Numerical methods for this are outlined in [13] and have been carried out for the linearized 
Navier-Stokes equations. They are, however, expensive and cannot be guaranteed to locate 
all critical points. It is, therefore, of interest to identify a priori problems for which s* = 0 
is a solution of (11-13). We find that sufficient conditions can be given if we make the 
following partial symmetrizability assumption. 

that u = Su 
(25) 

D>0 , M > O, Doj=0, j = 1 , . . . n — 1 . (26) 

Furthermore , the boundary conditions on dQ either take the form: 

- du 

T dn + “ = °’ < 27 ) 

or 

du ~ 

Tn + Ti = °- < 28 > 

with T symmetric and positive semide finite. 


Assumption 1 There exists a smooth matrix S(y) with smooth inverse such 
satisfies: 

~ du ~ du - _ - 

M T A — — T Cu — 72 V ■ Z) V u , 
dt dx 

where the matrices M , A y C and D are symmetric and: 


We note that Assumption 1 will hold for a number of interesting physical problems in- 
cluding some advection-diffusion equations as well as linearized advection-reaction-diffusion 
equations with parallel base flows. It does not, however, hold for the Navier-Stokes equa- 
tions linearized about a parallel flow. In fact s* ^ 0 may be important in our asymptotic 
analysis of linearized viscous flow both for very low and very high Reynolds numbers [12]. 

We now turn to the analysis of (5). In what follows we’ll drop the tildes and assume 
that (1- 4) has been given in the partially symmetrized form above. If we can show that 
A/(s) is real for real s in an interval containing zero, than we will have shown that the real 
5-axis is locally a curve of real group velocities. Note that the eigenvalue problem is not 
self-adjoint due to the dependence of Q on A. Nonetheless, by a simple adaptation of an 
argument due to Berestycki and Nirenberg [3], we are able to prove: 
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Theorem 1 Suppose Assumption 1 holds and, if the boundary condition is given by (28), 
the matrix f dQ T + f n C is positive definite. Then there exists 6 < 0 such that solutions, 
A(s), of (5) are real when s is real and s > 6. 


Proof: Following [3] we write A = p + w and by an integration by parts obtain: 

B(v, A) + s f v h Mv = XQ(v,p), (29) 

Jn 

. \ ( v h Tv ' 

wvm + ^ +j L w 

\ ' dn dn > 

Q(v,p)= f v H (2pVDoo — A)v. ( 31 ) 

J Q 

As Q and B are real, if 5 is real A can have a nonzero imaginary part only if both sides of 
(29) are zero. The hypothesis of the theorem, moreover, implies the existence of 6 < 0 such 
that for s > 8 the left-hand side is positive, guaranteeing that A is real and of the same sign 

Equation (29) can further be used to establish a connection between the sign of A/(0) 
and both its large s behavior and the direction of propagation of the corresponding wave 
group: 

Theorem 2 Suppose the conditions of Theorem 1 hold and M > 0. Then l G A f and 
A{(0) < 0 if and only if A/(0) < 0. 


Proof: Consider any nonsingular solution branch A/(s) as s varies along the nonnegative 
real axis. As the left-hand side of (29) is strictly positive, the right-hand side can never be 
zero, so A; must be of one sign. Moreover, if A/ < 0 then Q < 0. Differentiating (5) with 
respect to s, multiplying by v T and integrating we obtain a general formula for AJ: 


/ Mv = A; / (2\iVDoo — A)v. (32) 

Jn Jn 

For A i and v real the factor multiplying A{ is simply Q so that the sign of A} is the same as 
that of Q and, hence, that of A;. 

Note that the conditions be relaxed to allow Aj(0) = 0 so long as £7(u(i/; 0), 0) ^ 0. Then 
the sign of Q determines the sign of A{ and A;. It should also be noted that these theorems, 
though establishing the contribution of the neighborhood of s = 0 to the asymptotics, do 
not rule out the existence of imaginary s contributions. That is, we have been unable to 
prove that real group velocities cannot occur for imaginary s even with the assumptions 
made above. 


3 Construction of Asymptotic Boundary Conditions 

Now suppose we are interested in the numerical solution of (1-4), or, more generally, of a 
problem for which these equations represent a far-ficld approximation. In order to restrict 
the computational domain to a finite region, an artificial boundary is typically introduced. 


say at x = r. In order to close the system, additional boundary conditions must be imposed. 
That is, we approximate u by u f which satisfies (1-4) for (x,y) 6 [0,r] x ft as well as: 

= o. (33) 

In the past fifteen years an extensive literature has developed discussing the problem of 
choosing the operator B , though mainly in the difficult case of hyperbolic problems. There 
are three basic criteria which it should ideally satisfy: 

1. The finite domain problem is well-posed. 

2. Bu is small in an appropriate norm. 

3. The efficient numerical solution of the finite domain problem is possible. 

The first two criteria are the specialization of the usual notions of stability and consistency 
to the problem of numerical boundary conditions, and can be used to derive error estimates 
[13]. Techniques for establishing the first criterion are well-known [7], [17], but the second 
is far more difficult to satisfy. In particular, exact conditions, i.e. those for which Bu = 0 
whenever u solves (1) and (4), are invariably nonlocal in both space and time and, hence, 
violate the third criterion. Asymptotic analysis can provide a means for constructing local, 
or at least local in time, operators which still satisfy the ‘consistency’ criterion. 

The Laplace transforms of ^ and qi are related by: 

- 7 ^{x,y;s) = \i{s)qi(x,y;s). (34) 

If the assumptions of the previous section are satisfied, the asymptotic analysis involves the 
restriction of the transforms to a neighborhood of s = 0. An asymptotic expansion of the 
x derivative may then be obtained by replacing A /($) by its Taylor series about 0: 

A/(s) A/(0) + Aj(0)s + . . . . (35) 

This leads to the local relationship: 

If ~ («“> + a :< 0 4) «■ w 

These may be substituted into (9) to finally obtain a condition on tq. The time derivative is 
brought outside the integral to further simplify the expression. This involves the neglect of 
terms from the limits of integration which should be exponentially small. The asymptotic 
boundary condition we propose is, then, given by: 

§7 = (mo) + A I(0)^) (37) 

The construction of B now hinges on the number of modes, tt/, which are important in 


the 

asymptotic development of u. In tl 

ie simplest case there is a single mode, ss 

IV 1=1, 

sati 

sfying: 




»c rA '<°>, 1/ 1. 

(38) 

The 

n we may take: 




0 

0 



H Ox 

» |cS 

o 

^r 

t 

o 

I 

(3!)) 


0 


In [13] it is shown that this procedure leads to an error estimate of the form: 

e T M°) . . 

||u - «;ll < — H“ll. < 40 > 

where K is a constant. In the most interesting cases Ai(0) is smalt so that the 7 decay, 
which is a direct result of the use of asymptotic conditions, is important. 

More generally we may assume that A/(0) < A/+i(0) and, for some J . 

e ^M°) 1(°) # (41) 


If J is small, for example 2 or 3, a Jth order operator B may be constructed as the product 
of first order conditions: 


*-n (^-Mo)-Ai(°)|)- (42) 

For J larger, a nonlocal in space condition is likely to be more useful. Let /C;, i — 1, • • • ,J 
be linear functionals on the appropriate space of functions on with the property: 

£,t>i(-;0) = 6u, i,/=l,...,J. ( 43 ) 


The existence of such functionals is guaranteed by the linear independence of the eigenfunc- 
tions vi(y;0). An asymptotic boundary operator is then given by: 


„ d . .8 


(44) 


where 

Ao« = EA t (0)(AC l u)u.( 3 /;0), ( 45 ) 

t — 1 

J 

Aiu = 53A;(0)(/Citx)i;i(y ; 0). ( 46 ) 

t=i 

Product conditions with J = 2 have been used in [11] and [12], but the nonlocal condition 
(44) is as yet untried. 


3.1 Behavior for T> Small 

The study of waves in dissipative systems with small dissipation, in our case for small P, 
is of particular interest. For the problem of boundary conditions it provides a potential 
link between our results and the hyperbolic theory. In a recent paper Halpern [15] has 
studied this question for incompletely parabolic perturbations of linear hyperbolic systems 
with constant coefficients. The conditions derived there involve viscous corrections to the 

hyperbolic conditions of Engquist and Majda [8]. 

Relaxing the the assumptions of the previous section to allow A(0) = 0, we seek solutions 
of (5) with bounded decay rates. Thus the term proportional to PA 2 can be ignored to 
leading order. We distinguish between two separate cases: 

(la, sc 1: There are no smooth functions v(y) satisfying (4) and Cv — 0. 

Case 2 : ( ' = (). . 

(We are excluding the more complicated case of C f 0 but violating the case 1 assumption.) 


I 


Under case 1 we must consider solutions of the singular perturbation problem: 

n_1 


V E + Cu = -XAV. 


( 47 ) 


This is rather difficult to analyze, and we won’t discuss it further. Note, however, that we 
expect the decay rates to be bounded away from zero independent of V. 

Under case 2, on the other hand, we find many eigenvalues which approach 0 with P. 
Introducing the scaling A = T > A we consider the nonsingular reduced problem: 


n — i 


- E 


_d_ 

dyi 


D; 


dv_ 

dyj 


-XAv. 


(48) 


All eigenvalues of the reduced problem which are negative are of importance in the asymp- 
totic boundary conditions for D sufficiently small, as their decay rates are approaching 
zero. This suggests the use of the nonlocal conditions of the preceding section. If A > 0, 
which we may associate with all characteristics pointing out of the computational domain, 
all eigenvalues are negative and the reduced equation itself may be used as a boundary 
condition: 


, , dti . du 
M ~di + a Tx 


n— 1 


v y —d — 


(49) 


If A < 0, corresponding to all characteristics pointing in, then there are no negative eigen- 
values bounded independent of V. Dirichlet conditions may then be imposed on u. The 
most interesting case occurs when A has eigenvalues of both sign. Then the boundary con- 
dition will involve a coupling between incoming and outgoing characteristics of the related 
hyperbolic problem. 

It should be noted that we do not expect the error estimates to be uniformly valid as 
T) approaches 0. For hyperbolic problems curves of real group velocity will coincide with 
the imaginary s axis, which indicates nonuniformity in the steepest descent calculations. 
That is, we cannot expect only the low frequencies to be present for moderate r. It would 
be of interest in this regard to extend our method to viscous perturbations of the wave 
equation and study its connections with the recent work of Engquist and Halpern [9]. They 
consider low frequency corrections to the Engquist-Majda conditions and study the long 
time behavior of the error. 


3.2 Applications to the Navier-Stokes Equations 

The Navier-Stokes equations provide a most important example of a dissipative system 
with propagating solutions and were the prime motivation for the work described here. The 
detailed application of our technique to the derivation of asymptotic boundary conditions 
for the incompressible equations linearized about plane parallel flow is described in [12]. 
I lie eigenvalue problem (5) is then given by the well-known Orr-Sommerfeld equation of 
hydrodynamic stability. As the equations do not satisfy Assumption 1, the results guar- 
anteeing s* rr 0 cannot be applied. Nonetheless, for large Reynolds number, we find that 
an analogue to (-18) can be constructed which does lead to slowly decaying propagating 
modes near s* - 0. A boundary condition based on two of these modes was utilized in the 
simulation of vortical disturbances to Poiseuille flow and was quite successful for Reynolds 
numbers in tin* hundreds. For higher Reynolds numbers a condition such as (-19) is likely to 


8 


be effective. For flows with a critical Reynolds number, however, solutions with imaginary 
s* will eventually become dominant. Our analysis suggests that the nonlocal condition (44) 
based on solutions of the full Orr-Sommerfeld equation may then be useful, but they have 
not yet been tried. 

Presently under consideration are applications to compressible flows and to nonisother- 
mal incompressible flows. Also of note are the detailed studies by Halpern and coworkers 
[16], [15] of boundary conditions for both the compressible and incompressible Navier-Stokes 
equations linearized about a uniform flow. 


4 Traveling Wave Solutions in Cylindrical Domains 

Among the most active areas in the study of wave propagation in dissipative systems has 
been the analysis of planar traveling wave solutions to reaction-diffusion equations. Appli- 
cations have been made to models of flame propagation, nerve conduction and the internal 
structure of shocks in viscous conservation laws to name a few. Simple examples include 
Fisher’s equation, which is the diffusion equation plus the nonlinearity F = u(l - u), where 
a continuum of waves is found connecting a linearly stable and linearly unstable homoge- 
neous state, and the bistable equation, F = u(u - a)(l — u), where a unique wave connects 
two stable states. 

A natural, nontrivial generalization of these solutions to higher space dimensions comes 
from the consideration of waves propagating in cylindrical domains. In general we consider. 

^ = T, S- D nP- + F(u,Vu,y), (x,y) 6 (- 00 , 00 ) x ft. (50) 

dt 0 dy t dyj 

Here again we identify j/o with x and suppose that neither D,y nor F depend explicitly 
on x. Boundary conditions (4) are also imposed. We further suppose that (at least) two 
solutions independent of x and t, u±(y) exist. A traveling wave solution of (50) takes the 
form u = w(x — ct,y) with: 

lim z ^± (X> w(z } y) = u±(y). ( ol ) 

It satisfies the elliptic equation: 


n— 1 




d _ dw dw . \ n 

Dij— + c— + F(w, Vin, y) = 0. 
dyj dz 


(52) 


Only a few studies of solutions to (52) have appeared in the literature. Berestycki and 
Nirenberg [3], [2] have established monotonicity results for general scalar problems. Gardner 
[10] proves the existence of a nonplanar wave for the bistable equation in a channel. In joint 
work with Buonincontri [5], we have numerically computed traveling wave solutions to the 
Frank- Kamenetski equations in a channel and also developed some stability results. As the 
waves connect stable and unstable states, this may be thought of as a multidimensional 
problem of Fisher type. 
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4.1 A Model of Flame Propagation in a Channel with Cold Walls 

Consider the following system of reaction-diffusion equations modeling combustion in the 
channel {x,y) 6 (- 00 , 00 ) X [0, 1]: 


dY 
’ dt 


dT , -e 

-gt-VT.vr* r. 

(53) 

. _£-iy 2 y = -VpYe^r, 

(54) 

pT = k, constant, 

(55) 

dY 

T 0 ), ~dn~ Q ' at y = 0,1. 

(56) 


dT „ 

Here T is the temperature, Y is the mass fraction of the reactant, C is the Lewis number and 
8 is the activation energy. The model has heat loss at the channel walls, which precludes 
the existence of planar traveling wave solutions. 

This thermo-diffusive model ignores convection by the fluid flow induced by thermal 
expansion. Most authors, in using such a model, take p constant rather than proportional 
to T . In order to define the traveling wave we must find unburnt and burnt limiting states 
at ± 00 . The burnt state solution has no reactant present and is at the ambient temperature: 


Y = 0, T = r o . 


(57) 


Note that in the study of plane flames the burnt boundary is hot, T > T 0 . The structure of 
the flame we compute will be much different, exhibiting a nonmonotonic temperature profile 
rising to a maximum in the reaction zone then slowly decaying due to the heat loss. There 
is no unburnt equilibrium due to the form we’ve chosen for the reaction term. However, 
supposing the dimensionless ratio of ambient temperature to activation energy to be small, 
c = ^ < 1, we can compute a slowly varying unburnt state: 

Y = Y o {tt) + cY l {it,y) + 0{c' i ), T = T 0 + eT 1 (€t,y) + 0(e 2 ). (58) 


Among the questions one would like to answer about the propagating flames are: 

1. For what values of the parameters do traveling waves exist? Is there extinction for suf- 
ficiently large heat loss? (Evidence for extinction will be nonexistence of propagating 
solutions.) 

2. How does the shape and speed of the waves vary with the parameters and compare 
with planar results? 

A detailed discussion of the derivation of the model equations can be found in the doctoral 
dissertation of Buonincontri [4], In a recent paper, Benkhaldoun, Larrouturou and Denct [1] 
present a numerical study of essentially the same problem using a slightly different thermo- 
diffusive system. Rather than directly solving the equations defining the traveling wave, 
they solve an initial-boundary value problem. 
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4.2 Numerical Methods and Preliminary Results 

In this section we shall describe our numerical methods for computing the traveling wave, 
which are of general interest, and also present some preliminary results. All the compu- 
tations are due to Buonincontri. A more detailed and definitive study will be published 
elsewhere when complete. 

The fundamental numerical difficulty is the proper imposition of boundary conditions 
at artificial boundaries, z = z±, introduced to limit the computational domain, and the 
determination of the unknown wave speed c. We are guided by the general theory o 
asymptotic boundary conditions for nonlinear elliptic boundary value problems given in 
[14]. For ordinary differential equations the numerical problem has also been considered by 
Doedel and Friedman [ 6 ]. They obtain impressive theoretical results, but the application of 
their technique to partial differential equations would be too expensive. 

An asymptotic analysis of the traveling wave as z -* ±oo yields: 


( T -T 0 - eTi 
^ y - y 0 - cTi 


/V 


7+e 


-A + z 


v+(y) 

™+(y) 


z — ► CO, 


(59) 


( T -To 

v y 


7_e 




v-(y) 

w-(y) 


z — > — 00 . 


(60) 


The exponential decay rates, ±A, and corresponding eigenfunctions are defined by the 
eigenvalue problem (5) for the equations linearized about the asymptotic states. It should 
be noted that they depend on the unknown wave speed c. From these expansions we derive 


the asymptotic boundary conditions: 



The traveling wave problem (52) clearly cannot possess unique solutions due to the trans- 
lation invariance of the wave. To get uniqueness, we must impose an additional phase 
condition which we do by fixing 7 +: 


j\v + {T - {T 0 + e7i)) + to+(y - (y o + <Y\))}dy = k. (63) 

Jo 

This additional equation balances the additional unknown, c. In [5], where a stable-unstable 
connection is computed, the use of the theory of [14] leads to one less boundary condition 
at z + corresponding to the fact that waves exist for a range of speeds. 

We approximate (52) using centered finite differences and solve for the wave and c 
simultaneously using Newton’s method. Note that the et dependence of the wave only 
comes in, at leading order, through the boundary condition. In Figure 4.2 we plot a typical 
profile. At present we have only investigated a small part of parameter space with C > 1. 
We have found a variety of interesting phenomena including extinction for a large enough 
as well as an apparent limit point as £ is varied. 
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